Прежде чем перейти к рассмотрению модели ARIMA, познакомимся сначала с двумя другими моделями: скользящего среднего и моделью авторегрессии.

Модель скользящего среднего MA($q$)

Модель скользящего среднего порядка $q$ или просто MA($q$) предполагает следующую зависимость данных:

$$y_t\ =\ \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q},$$

где $y_t$ — стационарный ряд со средним $\mu$, а $\varepsilon_t$ — гауссовский белый шум, то есть $\varepsilon_t \sim \mathcal{N}(0, \sigma^2)$ и независимы.

По сути наш ряд $y_t$ выражается через сумму некоторого фиксированного среднего $\mu$, значения белого шума в текущий момент времени $\varepsilon_t$ и не более $q$ предыдущих значений белого шума, домноженных на некоторые коэффициенты, которые являются параметрами модели.

Рассмотрим некоторые свойства модели MA($q$). Как уже было упомянуто выше, ряд $y_t$ будет являтьcя стационарным со средним $\mu$. Найдем также $\mathsf{D}y_t$. Воспользовавшись свойством независимости для $\varepsilon_t$, можем заключить, что

$$\mathsf{D}y_t = \left(1 + \theta_1^2 + \dots + \theta_q^2\right)\sigma^2.$$

Посчитаем автоковариационную функцию для ряда $y_t$, то есть найдем значение $cov(y_t, y_{t + \tau})$. Легко понять, что если $\tau > q$, то $cov(y_t, y_{t+\tau})$ = 0, т.к. $\varepsilon_t$ независимы. Если же $\tau \leq q$, то тогда

$$cov(y_t, y_{t+\tau}) = \left(\theta_{\tau} + \theta_1\theta_{\tau + 1} + \dots + \theta_{q-\tau}\theta_{q}\right)\sigma^2.$$

Записав более компактно, можем получить:

$$cov(y_t, y_{t+\tau}) = \begin{cases} \sigma^2\sum_{j=0}^{q-\tau}\theta_j\theta_{\tau + j}, & \tau \leq q; \\ 0 & \tau > q; \end{cases}$$

где $\theta_0 = 1$. Из посчитанных значений для дисперсии и ковариационной функции, можете попробовать получить выражение и для автокорреляционной функции. Ее особенностью будет как раз равенство нулю на лаге, превосходящим $q$.

Посмотрим на визуализацию:

Модель авторегрессии AR($p$)

Модель авторегрессии для временного ряда можно записать следующим образом:

$$y_t\ =\ \alpha + \varphi_1 y_{t-1} + ... + \varphi_p y_{t-p} + \varepsilon_t,$$

где $y_t$ — стационарный ряд, а $\varepsilon_t$ — гауссовский белый шум, то есть $\varepsilon_t \sim \mathcal{N}(0, \sigma^2)$ и независимы. Отметим, что, вообще говоря, для стационарности нужны некоторые условия на коэффициенты $\varphi_1, ..., \varphi_p$.

По сути наш ряд $y_t$ выражается через сумму некоторого фиксированного числа $\alpha$, значения белого шума в текущий момент времени $\varepsilon_t$ и не более $p$ предыдущих значений этого же ряда, домноженных на некоторые коэффициенты, которые являются параметрами модели.

Другими словами, модель AR($p$) — это модель $\textit{линейной регрессией}$ для которой

• Таргет: $y_t$ — значение ряда в момент времени $t$
• Признаки: $y_{t-1}, ..., y_{t-p}$ — значения ряда в предыдущие моменты времени

Введем $L$ — оператор сдвига, обладающий следующими свойствами:

• применение $L$ к ряду дает предыдущее значение этого же ряда: $Ly_t = y_{t-1}.$
• применение $L$ к белому шуму дает предыдущее значение шума: $L\varepsilon_t = \varepsilon_{t-1}.$
• применение $L$ к константе — это константа: $Lc = c.$

Оператор $L$ иногда называют также лаговым оператором. Можно рассматривать функции от оператора сдвига, например, кратное применение оператора $L$: $L^2 y_t = L(L y_t) = L(y_{t-1}) = y_{t-2}$ или $L^{-1} y_t= y_{t+1}$. Для записей некоторых моделей временных рядов будет удобно использовать лаговый многочлен:

$$\varphi(L) = \sum_{i=1}^p \varphi_i L^{i}$$

Обратным к оператору $\varphi(L)$ называют оператор $\varphi^{-1}(L)$ такой, что:

$$\varphi(L)\varphi^{-1}(L)y_t = \varphi^{-1}(L)\varphi(L)y_t=y_t$$

Так, например, для $\lvert \varphi \rvert < 1$ можно заключить, что:

$$\frac{1}{1 - \varphi L} = \left(1-\varphi L \right)^{-1} = \sum_{i=1}^{\infty}\varphi^{i}L^{i}$$

Рассмотрим модель AR($p$):

$$y_t\ =\ \alpha + \varphi_1 y_{t-1} + ... + \varphi_p y_{t-p} + \varepsilon_t$$

С помощью оператора сдвига ее можно представить в следующем виде:

$$a(L) y_t\ =\ \alpha + \varepsilon_t,$$

где $a(z) = 1 - \varphi_1 z - ... - \varphi_p z^p$ — характеристический полином.

Сформулируем пару важных утверждений:

• Любой стационарный (в широком смысле) процесс представим в виде $MA(\infty)$, то есть в виде модели скользящего среднего с неограниченным количеством слагаемых (конечное или бесконечное число). Этот результат так же известен как теорема Волда о декомпозиции временного ряда.
• Модель $AR(p)$ задает стационарный временной ряд $\Longleftrightarrow$ все комплексные корни $a(z)=0$ лежат вне единичного круга.

Приведем пояснение второго утверждения. В самом деле, пусть $z_1, ..., z_p$ — все его комплексные корни (их ровно $p$ с учетом кратности), тогда справедливо представление:

$$a(z) = (z-z_1)...(z-z_p) = z_1...z_p \left(1 - \frac{z}{z_1}\right) ... \left(1 - \frac{z}{z_p}\right)$$

Тогда при представлении временного ряда в виде

$$y_t\ =\ \frac{\alpha + \varepsilon_t}{a(L)}$$

и дальнейшего его разложения на простые дроби возникнут слагаемые вида

$$\frac{\varepsilon_t}{1 - \frac{L}{z_j}}.$$

Если при этом $z_j$ лежит внутри единичного круга или на его границе, то соответствующий ряд будет расходящимся. На самом деле, случай $z_j=1$ мы в дальнейшем учтем.

В качестве примера рассмотрим подробнее модель $AR(1)$.
Зависимость имеет вид $y_t = \alpha + \varphi y_{t - 1} + \varepsilon_t$, где $\varepsilon_t \sim \mathcal{N}(0, \sigma^2)$. Для данного ряда можно выписать следующие свойства:

• Уравнение $1 - \varphi z = 0$, имеет корень $\lambda = 1/\varphi$.
• Тем самым, $AR(1)$ стационарен $\Longleftrightarrow$ $\lvert\varphi\rvert < 1$. Кроме того, чем меньше $\varphi$, тем предыдущее значение ряда вносит меньший вклад в текущее значение.
• Если ряд стационарен, то:
  • $\mathsf{E} y_t = \frac{\alpha}{1-\varphi}$
  • $\mathsf{D} y_t = \dfrac{\sigma^2}{1-\varphi^2}$
  • $cov (y_t, y_{t - h}) = \varphi^h \cdot \dfrac{\sigma^2}{1-\varphi^2}$.

Разберем первое равенство, остальные получаются аналогично. Возьмем математическое ожидание в уравнении ряда

$$\mathsf{E} y_t = \alpha + \varphi \mathsf{E} y_{t - 1} +\mathsf{E} \varepsilon_t$$

Поскольку ряд стационарен, то его математическое ожидание не меняется во времени, а для белого шума математическое ожидание равно нулю. Тем самым мы получаем уравнение на $m = \mathsf{E} y_t$, откуда следует доказываемая формула.

Таким образом, в зависимости от значения $\varphi$ мы можем получить следующие результаты:

• Если $\lvert\varphi\rvert < 1$, то $y_t = \mu + \sum\limits_{j = 0}^{\infty}\varphi^j\varepsilon_{t - j}$ — представление ряда в виде MA($\infty$).
• Если $\lvert\varphi\rvert = 1$, то $AR(1)$ — это случайное блуждание.
• Если $\lvert\varphi\rvert > 1$, то $AR(1)$ — экспоненциально растущий процесс.

Посмотрим на визуализацию.

• В первом случае мы имеем модель $y_t = - 0.5 y_{t - 1} + \varepsilon_t$, отрицательный коэффициент является следствием больших колебаний ряда.

• Во втором случае модель $y_t = 0.9 y_{t - 1} + \varepsilon_t$, большой положительный коэффициент делает ряд менее шумным.

• В третьем случае показано несколько рядов вида случайного блуждания $y_t = y_{t - 1} + \varepsilon_t$, что соответствует случаю $\varphi=1$.

• В четвертом случае показан экспоненциальный процесс $y_t = 1.1 y_{t - 1} + \varepsilon_t$, на графике шум уже не заметен из-за масштаба.

На немного вернемся к модели MA($q$). Чуть выше мы выяснили, что при некоторых условиях на коэффиценты $\varphi$ временной ряд модели AR($p$) будет стационарным, а значит имеет представление в виде MA($\infty$). На самом деле, модель скользящего среднего порядка $q$ тоже можно представить с помощью оператора $L$ следующим образом:

$$y_t\ =\ \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q}\ =\ \mu + b(L)\varepsilon_t$$

где $b(z) = 1 + \theta_1 z + ... + \theta_q z^q$ — характеристический многочлен. Для простоты изложения пусть $\mu = 0$. Важным при такой записи оказывается понятие обратимости, то есть представления в виде

$$\varepsilon_t = b^{-1}(L)y_t,$$

которое означает, что ряд можно представить в виде бесконечной авторегрессионной модели.
Здесь, как и в рассуждениях выше, можно заключить, что временной ряд $y_t$ обратим, если все комплексные корни $b(z) = 0$ лежат вне единичного круга.

Модель ARMA($p, q$)

Модель ARMA($p, q$) по сути является суммой моделей $AR(p)$ и $MA(q)$, иначе говоря, модель есть сумма нескольких предыдущих значений ряда и нескольких предыдущих значений белого шума с некоторым коэффициентами.

$$y_t\ =\ \alpha + \varphi_1 y_{t-1} + ... + \varphi_p y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q}$$

Эквивалентную запись ряда в терминах оператора сдвига можно получить, рассмотрев два многочлена

$$a(L) y_t = \alpha + b(L) \varepsilon_t$$

или

$$y_t = \mu + \frac{b(L)}{a(L)} \varepsilon_t,$$

где $a(z) = 1 - \varphi_1 z - ... - \varphi_p z^p,$ и $b(z) = 1 + \theta_1 z + ... + \theta_q z^q$.
Заметим, что во втором представлении константа $\alpha$ заменена на $\mu = \mathsf{E} y_t$. На самом деле, стационарность такого ряда будет определяться только его AR($p$) компонентой, то есть значениями коэффициентов $\varphi$, так ряд в модели MA($q$) всегда является стационарным.

Модель ARIMA($p, d, q$)

Модель ARIMA($p, d, q$) — это расширение моделей типа ARMA на нестационарные временные ряды, которые однако могут стать стационарным после применениея процедуры дифференцирования ряда. Модель ARIMA($p, d, q$) для ряда $y_t$ определяется как модель ARMA($p, q$) для ряда разностей порядка $d$ ряда $y_t$.

• Разность порядка 1: $y_t - y_{t-1} = (1 - L) y_t$.
• Разность порядка 2: $(1 - L)^2 y_t = (1 - L) (y_t - y_{t-1}) = (y_t - y_{t-1}) - (y_{t-1} - y_{t-2}) = y_t - 2y_{t-1} + y_{t-2}$.

Получаем формулу модели ARIMA:

$$a(L) (1 - L)^d y_t = \alpha + b(L) \varepsilon_t$$

или

$$(1 - L)^d y_t = \mu + \frac{b(L)}{a(L)} \varepsilon_t.$$

То есть многочлен $\widetilde{a}(z) = a(z) (1-z)^d$ имеет $d$ единичных корней.
Тем самым такая модель позволяет учесть нестационарности, в частности, тренд.

В качестве примера рассмотрим процесс случайного блуждания:

$$y_t = y_{t-1} + \varepsilon_t,$$

где $\varepsilon_t$ — белый шум. Как уже упомяналось ранее, такой ряд не является стационарным. Однако, если мы применим операцию дифференцирования, то можем перейти к новому, уже стационарном ряду $y'_t = y_t - y_{t-1}$, который можно записать в виде:

$$y'_t = \varepsilon_t$$

Частичная автокорреляция

Для модели скользящего среднего порядка $q$ мы выяснили, что значения автокорреляционной функции для такого ряда оказывается равной нулю после лага $q$. Эта особенность позволяет использовать автокорреляционную функцию для определения порядка модели скользящего среднего. Возникает разумный вопрос, как оценить порядок $p$ для модели AR($p$)? Здесь оказывается полезным понятие частичной (частной) автокорреляционной функции.

Частичная автокорреляция (PACF) — корреляция ряда с собой после снятия линеной зависимости от промежуточных значений ряда. Иначе говоря, мы хотим как-то учесть опосредованного влияние промежуточных значений ряда и оценить непосредственное влияние $y_{t - \tau}$ на $y_t$. Чуть более формально частичную автокорреляцию можно записать следующим образом:

$$\gamma_{\tau} = \begin{cases} corr(y_{t+1}, y_t), & \tau=1; \\ corr\left(y_{t+\tau} - y_{t+\tau}^{h-1}, y_t - y_t^{\tau-1}\right), & \tau\geqslant2, \end{cases}$$

где $y_t^{\tau-1}$ — линейная регрессия на $y_{t-1}, y_{t-2}, ..., y_{t - (\tau-1)}$:

• $y_t^{\tau-1} = \varphi_1 y_{t-1} + \varphi_2 y_{t-2} + ... + \varphi_{\tau-1} y_{t - (\tau-1)}$
• $y_{t+\tau}^{\tau-1} = \varphi_1 y_{t+\tau-1} + \varphi_2 y_{t+\tau-2} + ... + \varphi_{\tau-1} y_{t+1}$

Пример для $\tau=2$:

$$\gamma_{2} = corr\left(y_{t+2} - \varphi_1 y_{t+1}, y_t - \varphi_1 y_{t-1}\right)$$

где $\varphi_1$ — МНК-оценка в модели $y_t = \varphi y_{t-1}$.

Можно показать, что значение частиной автокорреляции для модели авторегресии AR($p$) будет ненулевой для лагов $\tau \leq p$ и равняться нулю для лагов $\tau > p$. Имеет место быть полная аналогия с автокорреляционной функцией и моделью MA($q$). Таким образом, исследование поведения автокорреляционной и частичной автокорреляционной функции может быть использовано для определения порядка $q$ модели скользящего среднего и порядка $p$ модели авторегрессии соответсвтенно.

Оценка коэффициентов в ARIMA

Пусть гиперпараметры $p, d, q$ фиксированы.&tab;В предположении, что $\varepsilon_t$ — гауссовский белый шум, в нашей модели мы можем выписать функцию правдоподобия $L_y(\theta, \varphi, \alpha) = p_{\theta, \varphi, \alpha}(y_1, ..., y_T),$ где $p_{\theta, \varphi, \alpha}(a_1, ..., a_T)$ — соместная плотность. Из-за того, что $\varepsilon_t$ имеют нормальное распределение, она будет иметь разумный вид. Соответственно, в качестве оценок параметров берется оценка максимального правдоподобия.

Для поиска начальных приближение для параметров $p$ и $q$ воспользуемся автокорреляционной и частичной автокорреляционной функцией.

• Начальное приближение $p$: последний значимый пик у PACF.
• Начальное приближение $q$: последний значимый пик у ACF.

Далее обычно используется поиск по сетке вокруг подобранных значений, минимизируя информационный критерий:

• $AIC = -2\ell^* + 2(p+q+1)$ — критерий Акаике;
• $AIC_c = -2\ell^* + \frac{2(p+q+1)(p+q+2)}{T-p-q-2}$ — критерий Акаике (короткие ряды);
• $BIC = -2\ell^* + (\log T - 2)(p+q+1)$ — Байесовский информационный критерий или критерий Шварца,

где $\ell^* = \ln L_y\left(\widehat{\theta}, \widehat{\varphi}, \widehat{\alpha}\right)$ — логарифм функции правдоподобия, $T$ — длина временного ряда.

Приведем некоторый план при применению модели ARIMA для прогнозирования временных рядов.

1. Анализ выбросов: замена нерелевантых выбросов на NA или усреднение по соседним элементам.

2. Стабилизация дисперсии (преобразования).

3. Дифференцирование, если ряд не стационарен.

4. Выбор пилотных $p$ и $q$ по PACF и ACF.

5. Вокруг этих параметров подбираем оптим. модель по $AIC$/$AIC_c$.

6. Пошаговое построение прогноза:
   — для $t \leqslant T$: $\varepsilon_t \Longrightarrow \widehat{\varepsilon}_t = y_t - \widehat{y}_t$;
   — для $t > T$: $\varepsilon_t \Longrightarrow 0$;
   — для $t > T$: $y_t \Longrightarrow \widehat{y}_t$.

7. Построение предсказательного интервала:

   — если остатки модели нормальны и гомоскедастичны (дисперсия постоянна), то строится теоретический предсказательный интервал

   $$\widehat \sigma^2(h) = \widehat \sigma^2 \left(1 + \sum\limits_{i = 1}^{h - 1} \widehat{\psi}_i^2\right)$$

   где $h$ — горизонт прогнозирования, $\widehat\sigma^2$ — оценка на дисперсию шума $\varepsilon_t$, $\widehat{\psi}_i$ — коэф. для ряда при его представлении в виде бесконечного процесса скользящего среднего. И $\widehat\sigma^2$, и $\widehat{\psi}_i$ могут быть выражены через оценки на параметры $\varphi$ и $\theta$.
   — иначе интервалы строятся с помощью бутстрепа.

Модели SARIMA и ARIMAX

Рассмотрим некоторые расширение модели ARIMA. Обобщение модели ARIMA на ряды с наличием сезонной составляющей назвается SARIMA. Пусть $s$ — известная сезонность ряда. Добавим в модель ARIMA($p, d, q$) компоненты, отвечающие за значения в предыдущие сезоны. Тогда модель SARIMA $(p, d, q)\times (P, D, Q)_s$ может быть записана следующим образом:

$$(1 - L)^d (1 - L^s)^D y_t = \mu + \frac{b(L) B(L^s)}{a(L) A(L^s)} \varepsilon_t,$$

где

$$a(z) = 1 - \varphi_1 z - ... - \varphi_p z^p,$$

$$b(z) = 1 + \theta_1 z + ... + \theta_q z^q,$$

$$A(z) = 1 - \varphi_1^s z - \dots - \varphi_P z^P,$$

$$B(z) = 1 + \theta^s_1 z + \dots + \theta_Q^s z^Q.$$

Параметр сезонного дифференцирования $D$, а также параметры $P, Q$ подбираются из тех же соображений, что и для $p, d, q$, но только с поправкой, что делается это с учетом сезонности $s$. ARIMAX — обобщение модели ARIMA, которая учитывает некоторые экзогенные факторы. Пусть $x_t \in \mathbb{R}^n$ — ряд регрессоров, известный до начала прогноза.

Простой вариант:

$$(1 - L)^d y_t = \mu + \sum_{i=1}^n \frac{\beta_i}{a(L)} x_t^i + \frac{b(L)}{a(L)} \varepsilon_t$$

Общий случай:

$$(1 - L)^d y_t = \mu + \sum_{i=1}^n \frac{u_i(L)}{v_i(L)} x_t^i + \frac{b(L)}{a(L)} \varepsilon_t$$

Пример: $x_t = I \{\text{в момент времени t праздник}\}$

Вышеуказанные модели можно объединить и получить SARIMAX $(p, d, q)\\times (P, D, Q)_s$:

$$(1 - L)^d (1 - L^s)^D y_t = \mu + \sum_{i=1}^n \frac{u_i(L)}{v_i(L)} x_t^i + \frac{b(L) B(L^s)}{a(L) A(L^s)} \varepsilon_t$$

Проведем аналогию с линейной регрессией. Это линейная по признакам модель, в которой

• Отклик: $y_t$ — значение ряда в моменты времени $t$
• Признаки:
  • $y_{t-1}, ..., y_{t-p}$ — значения ряда в предыдущие моменты времени
  • Значение ряда за предыдущие сезоны
  • Значения признаков в предыдущие моменты времени
  • Значения признаков в предыдущие сезоны
• Ошибка: сумма шума за предыдущие моменты времени и предыдущие сезоны.